In this script we will provide some intuition for harmonic regression, and walk through how to reconstruct temporally dynamic coefficients when fitting step selection (or other) models with harmonic terms. We will use the example of a model fitted with two pairs of harmonics. We will then use the coefficients to reconstruct the selection surfaces when we interact linear and quadratic terms with the harmonics.
Load packages
options(scipen=999)
library(tidyverse)
packages <- c("lubridate", "tictoc", "TwoStepCLogit", "beepr", "scales")
walk(packages, require, character.only = T)
Harmonics are sine and cosine functions are periodic functions that repeat once per cycle, \(T\). When combined, they can be used to model complex periodic patterns, and with many harmonics can create very flexible functions.
Harmonics can be considered similarly to splines or basis functions in that they are additive - if we add them all together we get a single, flexible function.
Let’s say we have a single pair of harmonics, \(s_1 = \sin(2\pi t / T)\), and \(c_1 = \cos(2\pi t / T)\), where \(t \in T\). For a yearly cycle for instance, \(t\) would be indexed as the day of the year, and \(T = 365\), although \(t\) does not need to be integer-valued, and can be arbitrarily fine. In the case of the model we are fitting here, \(t\) is the hour of the day, for \(T\) is the number of hours in a day, 24.
# create a sequence of hours
t <- seq(0,23,1)
T <- 24
# create the harmonic terms
sin1 <- sin(2 * pi * t / T)
cos1 <- cos(2 * pi * t / T)
# plot the harmonic terms
plot(t, sin1, type = "l", col = "red",
xlab = "Hour of the day", ylab = "Value",
main = "Harmonic terms", ylim = c(-2.5,2.5))
lines(t, cos1, col = "blue")
lines(t, rep(0,T), lty = "dashed")
Now what we can do is add them together to create a single function.
# create a single function
f_harmonic <- sin1 + cos1
# plot the function
plot(t, f_harmonic, type = "l", col = "black",
xlab = "Hour of the day", ylab = "Value",
main = "Harmonic function",
ylim = c(-2.5,2.5))
lines(t, rep(0,T), lty = "dashed")
It looks similar, although now the amplitude is greater (up to about 1.4 now), and that the peak has shifted to be an average of the sine and cosine functions (peak is now at 3 ( = 0 + 6 / 2). We can shift the location of the peak by giving one of the terms more influence.
# create a single function
f_harmonic <- 2*sin1 + cos1
# plot the function
plot(t, f_harmonic, type = "l", col = "black",
xlab = "Hour of the day", ylab = "Value",
main = "Harmonic function", ylim = c(-2.5,2.5))
lines(t, rep(0,T), lty = "dashed")
If we give more weight to the sine function, the peak will shift to the right, which is closer to the peak of the sine function, and vice versa. We can see it has also changed the amplitude, which now peaks above y = 2.
We can also add a constant term, which will just shift the entire function up or down.
# create a single function
f_harmonic <- 1 + 2*sin1 + cos1
# plot the function
plot(t, f_harmonic, type = "l", col = "black",
xlab = "Hour of the day", ylab = "Value",
main = "Harmonic function", ylim = c(-2.5,2.5))
lines(t, rep(0,T), lty = "dashed")
We can see that adding a constant, and by weighting the sin and cosine terms we can start to create a flexible function. If we only use 1 pair of harmonics we can only have one period, but if we use more pairs we can have more periods, and therefore more flexibility. Let’s add some more pairs of harmonics.
# create two more harmonic terms
sin2 <- sin(4 * pi * t / T)
cos2 <- cos(4 * pi * t / T)
# plot the harmonic terms
plot(t, sin2, type = "l", col = "red",
xlab = "Hour of the day", ylab = "Value",
main = "Harmonic terms", ylim = c(-2.5,2.5))
lines(t, cos2, col = "blue")
lines(t, rep(0,T), lty = "dashed")
Again we can add these together, and weight to shift the peak and change the amplitude.
# create a single function
f_harmonic <- 0.5*sin2 + 1.5*cos2
# plot the function
plot(t, f_harmonic, type = "l", col = "black",
xlab = "Hour of the day", ylab = "Value",
main = "Harmonic function", ylim = c(-2.5,2.5))
lines(t, rep(0,T), lty = "dashed")
Now we can start to add all of our components together to create flexible functions (and there is really no limit to this - infinitely flexible functions can be created with infinitely many harmonic terms, which is called a Fourier series - for some very cool intuition about Fourier series check out: https://www.youtube.com/watch?v=r6sGWTCMz2k).
# create a single function
f_harmonic <- 1 + 0.5*sin1 + 2*cos1 + 0.75*sin2 + 1.25*cos2
# plot the function
plot(t, f_harmonic, type = "l", col = "black",
xlab = "Hour of the day", ylab = "Value",
main = "Harmonic function")
lines(t, rep(0,T), lty = "dashed")
If you play around with the weights you can estimate some really funky functions, and the function will just get more flexible the more harmonic terms that you add.
# create more harmonic terms...
sin3 <- sin(6 * pi * t / T)
cos3 <- cos(6 * pi * t / T)
sin4 <- sin(8 * pi * t / T)
cos4 <- cos(8 * pi * t / T)
# create a single function
f_harmonic <- 0.5 + 0.5*sin1 + 1.5*cos1 + 2*sin2 + 1.25*cos2 + 0.25*sin3 + 0.5*cos3 + 0.1*sin4 + 2*cos4
# plot the function
plot(t, f_harmonic, type = "l", col = "black",
xlab = "Hour of the day", ylab = "Value",
main = "Harmonic function")
lines(t, rep(0,T), lty = "dashed")
Hopefully now you have some intuition about what we are doing with the harmonic terms in our step selection function, in that we are estimating the weights of the sine and cosine functions, which are the coefficients from our fitted models. Pretty cool!
Let’s have a look at some real data now, and fit a model with the harmonics.
Importing buffalo data for model fitting, and creating harmonic terms using the time at the end of the step.
buffalo_data_all <- read_csv(
"outputs/buffalo_popn_GvM_covs_ST_KDEmem1000_allOPTIM_10rs_2024-02-05.csv")
## Rows: 1082697 Columns: 26
## ── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (22): id, burst_, x1_, x2_, y1_, y2_, sl_, ta_, dt_, hour_t2, step_id_, y, ndvi_temporal, veg_woody, veg_herby...
## lgl (1): case_
## dttm (3): t1_, t2_, t2_rounded
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
buffalo_data_all <- buffalo_data_all %>%
mutate(t1_ = lubridate::with_tz(buffalo_data_all$t1_, tzone = "Australia/Darwin"),
t2_ = lubridate::with_tz(buffalo_data_all$t2_, tzone = "Australia/Darwin"))
buffalo_data_all <- buffalo_data_all %>%
mutate(id_num = as.numeric(factor(id)),
step_id = step_id_,
x1 = x1_, x2 = x2_,
y1 = y1_, y2 = y2_,
t1 = t1_,
t1_rounded = round_date(buffalo_data_all$t1_, "hour"),
hour_t1 = hour(t1_rounded),
t2 = t2_,
t2_rounded = round_date(buffalo_data_all$t2_, "hour"),
hour_t2 = hour(t2_rounded),
hour = hour_t2,
yday = yday(t1_),
year = year(t1_),
month = month(t1_),
sl = sl_,
log_sl = log(sl_),
ta = ta_,
cos_ta = cos(ta_),
# scale canopy cover from 0 to 1
canopy_01 = canopy_cover/100,
# recover the natural scale of previous space use density
kde_memory_density = exp(kde_memory_density_log),
# here we create the harmonic terms for the hour of the day
# for seasonal effects, change hour to day of the year (which is t),
# and 24 to 365 (which is T)
hour_s1 = sin(2*pi*hour/24),
hour_s2 = sin(4*pi*hour/24),
hour_s3 = sin(6*pi*hour/24),
hour_c1 = cos(2*pi*hour/24),
hour_c2 = cos(4*pi*hour/24),
hour_c3 = cos(6*pi*hour/24))
# to just select a single year of data
buffalo_data_all <- buffalo_data_all %>% filter(t1 < "2019-07-25 09:32:42 ACST")
buffalo_ids <- unique(buffalo_data_all$id)
First we create a data matrix to be provided to the model, and then we scale and centre the full data matrix, with respect to each of the columns. That means that all variables are scaled and centred after the data has been split into wet and dry season data, and also after creating the quadratic and harmonic terms (when using them).
This is where we interact the harmonic terms with our covariates. As we have already created the harmonic terms for the hour of the day (s1, c1, s2, etc), we just interact (multiply) these with each of the covariates, including the quadratic terms and movement parameters. As before, we create the data matrix with all quadratic and harmonic terms, and then scale the matrix by each column, and store the scaling and centering variables to reconstruct the natural scale coefficients (i.e. if we were to fit the model without centering and scaling).
months_wet <- c(1:4, 11:12)
buffalo_ids <- unique(buffalo_data_all$id)
# buffalo_data <- buffalo_data_all %>% filter(month %in% months_wet) # wet season
buffalo_data <- buffalo_data_all %>% filter(!month %in% months_wet) # dry season
buffalo_data_matrix_unscaled <- buffalo_data %>% transmute(
ndvi = ndvi_temporal,
ndvi_s1 = ndvi_temporal * hour_s1,
ndvi_s2 = ndvi_temporal * hour_s2,
ndvi_c1 = ndvi_temporal * hour_c1,
ndvi_c2 = ndvi_temporal * hour_c2,
ndvi_sq = ndvi_temporal ^ 2,
ndvi_sq_s1 = (ndvi_temporal ^ 2) * hour_s1,
ndvi_sq_s2 = (ndvi_temporal ^ 2) * hour_s2,
ndvi_sq_c1 = (ndvi_temporal ^ 2) * hour_c1,
ndvi_sq_c2 = (ndvi_temporal ^ 2) * hour_c2,
canopy = canopy_01,
canopy_s1 = canopy_01 * hour_s1,
canopy_s2 = canopy_01 * hour_s2,
canopy_c1 = canopy_01 * hour_c1,
canopy_c2 = canopy_01 * hour_c2,
canopy_sq = canopy_01 ^ 2,
canopy_sq_s1 = (canopy_01 ^ 2) * hour_s1,
canopy_sq_s2 = (canopy_01 ^ 2) * hour_s2,
canopy_sq_c1 = (canopy_01 ^ 2) * hour_c1,
canopy_sq_c2 = (canopy_01 ^ 2) * hour_c2,
slope = slope,
slope_s1 = slope * hour_s1,
slope_s2 = slope * hour_s2,
slope_c1 = slope * hour_c1,
slope_c2 = slope * hour_c2,
herby = veg_herby,
herby_s1 = veg_herby * hour_s1,
herby_s2 = veg_herby * hour_s2,
herby_c1 = veg_herby * hour_c1,
herby_c2 = veg_herby * hour_c2,
spatial_memory = kde_memory_density,
spatial_memory_s1 = kde_memory_density * hour_s1,
spatial_memory_s2 = kde_memory_density * hour_s2,
spatial_memory_c1 = kde_memory_density * hour_c1,
spatial_memory_c2 = kde_memory_density * hour_c2,
spatial_memory_sq = kde_memory_density ^ 2,
spatial_memory_sq_s1 = (kde_memory_density ^ 2) * hour_s1,
spatial_memory_sq_s2 = (kde_memory_density ^ 2) * hour_s2,
spatial_memory_sq_c1 = (kde_memory_density ^ 2) * hour_c1,
spatial_memory_sq_c2 = (kde_memory_density ^ 2) * hour_c2,
step_l = sl,
step_l_s1 = sl * hour_s1,
step_l_s2 = sl * hour_s2,
step_l_c1 = sl * hour_c1,
step_l_c2 = sl * hour_c2,
log_step_l = log_sl,
log_step_l_s1 = log_sl * hour_s1,
log_step_l_s2 = log_sl * hour_s2,
log_step_l_c1 = log_sl * hour_c1,
log_step_l_c2 = log_sl * hour_c2,
cos_turn_a = cos_ta,
cos_turn_a_s1 = cos_ta * hour_s1,
cos_turn_a_s2 = cos_ta * hour_s2,
cos_turn_a_c1 = cos_ta * hour_c1,
cos_turn_a_c2 = cos_ta * hour_c2)
buffalo_data_matrix_scaled <- scale(buffalo_data_matrix_unscaled)
# store the scaling and centering attribues to recover the natural scale coefficients
mean_vals <- attr(buffalo_data_matrix_scaled, "scaled:center")
sd_vals <- attr(buffalo_data_matrix_scaled, "scaled:scale")
scaling_attributes <- data.frame(variable = names(buffalo_data_matrix_unscaled),
mean = mean_vals, sd = sd_vals)
buffalo_data_scaled_2p <- data.frame(id = buffalo_data$id,
step_id = buffalo_data$step_id,
y = buffalo_data$y,
buffalo_data_matrix_scaled)
We add in all the terms in the model formula. We can see that we are adding all of these together, which allows us to add them together later to reconstruct the harmonic function.
formula_twostep <- y ~
ndvi +
ndvi_s1 +
ndvi_s2 +
ndvi_c1 +
ndvi_c2 +
ndvi_sq +
ndvi_sq_s1 +
ndvi_sq_s2 +
ndvi_sq_c1 +
ndvi_sq_c2 +
canopy +
canopy_s1 +
canopy_s2 +
canopy_c1 +
canopy_c2 +
canopy_sq +
canopy_sq_s1 +
canopy_sq_s2 +
canopy_sq_c1 +
canopy_sq_c2 +
slope +
slope_s1 +
slope_s2 +
slope_c1 +
slope_c2 +
herby +
herby_s1 +
herby_s2 +
herby_c1 +
herby_c2 +
spatial_memory +
spatial_memory_s1 +
spatial_memory_s2 +
spatial_memory_c1 +
spatial_memory_c2 +
spatial_memory_sq +
spatial_memory_sq_s1 +
spatial_memory_sq_s2 +
spatial_memory_sq_c1 +
spatial_memory_sq_c2 +
step_l +
step_l_s1 +
step_l_s2 +
step_l_c1 +
step_l_c2 +
log_step_l +
log_step_l_s1 +
log_step_l_s2 +
log_step_l_c1 +
log_step_l_c2 +
cos_turn_a +
cos_turn_a_s1 +
cos_turn_a_s2 +
cos_turn_a_c1 +
cos_turn_a_c2 +
strata(step_id) +
cluster(id)
As we have already fitted the model, we will load it here, but if the model_fit object file doesn’t exist, it will run the model fitting code. Be careful here that if you change the model formula, you will need to delete or rename the model_fit file to re-run the model fitting code, otherwise it will just load the previous model.
if(file.exists("outputs/model_twostep_2p_harms_dry.rds")) {
model_twostep_2p_harms <- readRDS("outputs/model_twostep_2p_harms_dry.rds")
} else {
tic()
model_twostep_2p_harms <- Ts.estim(formula = formula_twostep,
data = buffalo_data_scaled_2p,
all.m.1 = TRUE,
D = "UN(1)",
itermax = 10000)
toc()
# save model object
saveRDS(model_twostep_2p_harms, file = "outputs/model_twostep_2p_harms_dry.rds")
beep(sound = 2)
}
It will output a massive covariance matrix, which can be ignored for
our purposes at the moment. We are interested in the coefficients, which
are stored in the beta attribute of the model object.
model_twostep_2p_harms$beta
## ndvi ndvi_s1 ndvi_s2 ndvi_c1 ndvi_c2
## 1.18745765 0.39442696 0.01121123 -0.95521148 1.25104461
## ndvi_sq ndvi_sq_s1 ndvi_sq_s2 ndvi_sq_c1 ndvi_sq_c2
## -1.33022674 -0.24830951 -0.01598593 0.12688921 -0.46560189
## canopy canopy_s1 canopy_s2 canopy_c1 canopy_c2
## 0.27443845 0.16425934 -0.14215576 -0.03190371 0.25099088
## canopy_sq canopy_sq_s1 canopy_sq_s2 canopy_sq_c1 canopy_sq_c2
## -0.34884092 -0.09769922 0.04173728 -0.16166660 -0.16650068
## slope slope_s1 slope_s2 slope_c1 slope_c2
## -0.16917086 -0.01981839 -0.05097721 -0.03595167 -0.05467170
## herby herby_s1 herby_s2 herby_c1 herby_c2
## -0.02374414 -0.02758919 0.01137859 0.13312922 -0.04304990
## spatial_memory spatial_memory_s1 spatial_memory_s2 spatial_memory_c1 spatial_memory_c2
## 2.29925142 -0.48174644 -0.69537479 -0.61650328 0.48955057
## spatial_memory_sq spatial_memory_sq_s1 spatial_memory_sq_s2 spatial_memory_sq_c1 spatial_memory_sq_c2
## -1.39487590 0.17180133 0.55104843 0.60512613 -0.46486791
## step_l step_l_s1 step_l_s2 step_l_c1 step_l_c2
## -0.34761318 -0.10567643 -0.61073611 -0.15348490 -0.32241309
## log_step_l log_step_l_s1 log_step_l_s2 log_step_l_c1 log_step_l_c2
## 0.23121452 -0.13842027 -0.13823993 -0.23188181 0.06039618
## cos_turn_a cos_turn_a_s1 cos_turn_a_s2 cos_turn_a_c1 cos_turn_a_c2
## -0.01128474 -0.03790357 -0.18570605 -0.07010761 -0.02699471
Create a dataframe of the coefficients with the scaling attributes, and return the coefficients to their natural scale by dividing by the scaling factor (standard deviation).
As we can see, we have a coefficient for each covariate by itself,
and then one for each of the harmonics. These are the ‘weights’ that we
played around with above, and we reconstruct them in exactly the same
way. We also have the coefficients for the quadratic terms and the
interactions with the harmonics, which I have denoted as
ndvi_sq for instance. We will come back to these when we
look at the selection surfaces.
# creating data frame of model coefficients
coefs_clr <- data.frame(coefs = names(model_twostep_2p_harms$beta),
value = model_twostep_2p_harms$beta)
coefs_clr$scale_sd <- scaling_attributes$sd
coefs_clr <- coefs_clr %>% mutate(value_nat = value / scale_sd)
head(coefs_clr, 10)
## coefs value scale_sd value_nat
## ndvi ndvi 1.18745765 0.1429054 8.30939941
## ndvi_s1 ndvi_s1 0.39442696 0.2371431 1.66324473
## ndvi_s2 ndvi_s2 0.01121123 0.2372765 0.04724962
## ndvi_c1 ndvi_c1 -0.95521148 0.2426054 -3.93730465
## ndvi_c2 ndvi_c2 1.25104461 0.2425239 5.15843931
## ndvi_sq ndvi_sq -1.33022674 0.1019530 -13.04744654
## ndvi_sq_s1 ndvi_sq_s1 -0.24830951 0.1065570 -2.33029697
## ndvi_sq_s2 ndvi_sq_s2 -0.01598593 0.1076953 -0.14843668
## ndvi_sq_c1 ndvi_sq_c1 0.12688921 0.1105880 1.14740494
## ndvi_sq_c2 ndvi_sq_c2 -0.46560189 0.1097234 -4.24341508
Let’s have a look at herbaceous vegetation as it didn’t have a quadratic term.
Firstly, as above, we need a sequence of values that covers a full period (or the period that we want to construct the function over, which can be more or less than 1 period). The sequence can be arbitrarily finely spaced. The smaller the increment the smoother the function will be for plotting. When simulating data from the temporally dynamic coefficients, however, one would use the increment that relates to the data collection and model fitting (i.e. one hour in this case).
Here we’ll use a finer resolution than above for smoother plotting
hour_seq <- seq(0,23.9,0.1)
Now we can reconstruct the harmonic function using the formula that we put into our model by interacting the harmonic terms with each of the covariates, for a single covariate, let’s say herbaceous vegetation, this would be written down as:
\[ f = \beta_{herby} + \beta_{herby\_s1} \sin\left(\frac{2\pi t}{24}\right) + \beta_{herby\_c1} \cos\left(\frac{2\pi t}{24}\right) + \beta_{herby\_s2} \sin\left(\frac{4\pi t}{24}\right) + \beta_{herby\_c2} \cos\left(\frac{4\pi t}{24}\right) \]
We can precompute the harmonic terms to make it a bit neater (as we did for the design matrix before fitting our model)
t <- hour_seq
T <- 24
sin1 <- sin(2 * pi * t / T)
cos1 <- cos(2 * pi * t / T)
sin2 <- sin(4 * pi * t / T)
cos2 <- cos(4 * pi * t / T)
Now we just pull out the relevant coefficients (which are scalars (single numbers), like the weights we used above), multiply them by the harmonic terms (which are vectors representing the function through time), and add them together.
herby_harmonic_function <- coefs_clr$value[which(coefs_clr$coefs == "herby")] +
coefs_clr$value[which(coefs_clr$coefs == "herby_s1")] * sin1 +
coefs_clr$value[which(coefs_clr$coefs == "herby_c1")] * cos1 +
coefs_clr$value[which(coefs_clr$coefs == "herby_s2")] * sin2 +
coefs_clr$value[which(coefs_clr$coefs == "herby_c2")] * cos2
# plot the function
plot(t, herby_harmonic_function, type = "l", col = "black",
xlab = "Hour of the day", ylab = "Coefficient value (scaled)",
main = "Herbaceous vegetation")
lines(t, rep(0,length(t)), lty = "dashed")
From these we can see that buffalo select for herbaceous vegetation in the early morning and late afternoon, but have a reasonably strong avoidance of it in the middle of the day (and therefore a selection for woody vegetation). This aligns with what we know about buffalo behaviour (and the climate in Northern Australia’s tropical savannas), in that they are likely to be seeking shelter during this time due to high temperatures and sun.
Now that we know how the coefficients of the harmonics can be combined to form a temporally dynamic function, we can use a shortcut, which is matrix multiplication.
First we create a matrix of the harmonics, which is just the sin and cos terms for each harmonic, and then we can multiply this by the coefficients to get the function. When we use two pairs of harmonics we will have 5 coefficients for each covariate (linear + 2 sin and 2 cos), so there will be 5 columns in the matrix.
For matrix multiplication, the number of columns in the first matrix must be equal to the number of rows in the second matrix. The result will then have the same number of rows as the first matrix and the same number of columns as the second matrix.
Or in other words, if we have a 24 x 5 matrix of harmonics and a 5 x 1 matrix of coefficients, we will get a 24 x 1 matrix of the function, which corresponds to our 24 hours of the day (and the middle numbers be the same).
# we'll return the increments back to 1 hour for this example,
# just so everything is a bit clearer
hour_seq <- seq(0,23,1)
hour_harmonics_matrix <- as.matrix(data.frame("linear_term" = rep(1, length(hour_seq)),
"hour_s1" = sin(2*pi*hour_seq/24),
"hour_s2" = sin(4*pi*hour_seq/24),
"hour_c1" = cos(2*pi*hour_seq/24),
"hour_c2" = cos(4*pi*hour_seq/24)))
# now have a 24 x 5 matrix
head(hour_harmonics_matrix)
## linear_term hour_s1 hour_s2 hour_c1 hour_c2
## [1,] 1 0.0000000 0.0000000 1.0000000 1.00000000000000000000000
## [2,] 1 0.2588190 0.5000000 0.9659258 0.86602540378443870761060
## [3,] 1 0.5000000 0.8660254 0.8660254 0.50000000000000011102230
## [4,] 1 0.7071068 1.0000000 0.7071068 0.00000000000000006123032
## [5,] 1 0.8660254 0.8660254 0.5000000 -0.49999999999999977795540
## [6,] 1 0.9659258 0.5000000 0.2588190 -0.86602540378443870761060
We can now pull out the coefficients for each covariate and multiply by the matrix to get our temporally dynamic function. I’ll do this bit by bit initially.
Starting with the scaled NDVI covariate
# this subsets the coefficients for the covariate of interest
# (grepl uses string matching, and we include ndvi but exclude the quadratic term)
ndvi_coefs <- coefs_clr %>% filter(grepl("ndvi", coefs) & !grepl("sq", coefs))
ndvi_coefs
## coefs value scale_sd value_nat
## ndvi ndvi 1.18745765 0.1429054 8.30939941
## ndvi_s1 ndvi_s1 0.39442696 0.2371431 1.66324473
## ndvi_s2 ndvi_s2 0.01121123 0.2372765 0.04724962
## ndvi_c1 ndvi_c1 -0.95521148 0.2426054 -3.93730465
## ndvi_c2 ndvi_c2 1.25104461 0.2425239 5.15843931
# pull out the coefficients - this comes out as a column vector,
# which we'll convert to a matrix
ndvi_coefs_matrix <- matrix(ndvi_coefs %>% pull(value))
ndvi_coefs_matrix
## [,1]
## [1,] 1.18745765
## [2,] 0.39442696
## [3,] 0.01121123
## [4,] -0.95521148
## [5,] 1.25104461
# now we can multiply the matrix by the coefficients to get the function
ndvi_harmonic_function <- hour_harmonics_matrix %*% ndvi_coefs_matrix
# ndvi_harmonic_function
# plot the function
plot(hour_seq, ndvi_harmonic_function, type = "l", col = "black",
xlab = "Hour of the day", ylab = "Coefficient value",
main = "NDVI - linear term")
lines(hour_seq, rep(0,length(hour_seq)), lty = "dashed")
Now we can just repeat the process for all of the covariates, and put it all into a dataframe.
# back to the finer resolution for smoother plotting
hour <- seq(0,23.9,0.1)
hour_harmonics_matrix <- as.matrix(data.frame("linear_term" = rep(1, length(hour)),
"hour_s1" = sin(2*pi*hour/24),
"hour_s2" = sin(4*pi*hour/24),
"hour_c1" = cos(2*pi*hour/24),
"hour_c2" = cos(4*pi*hour/24)))
harmonics_df_2p <- data.frame(
"hour" = hour,
"ndvi" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("ndvi", coefs) & !grepl("sq", coefs)) %>%
pull(value)),
"ndvi_2" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("ndvi_sq", coefs)) %>% pull(value)),
"canopy" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("canopy", coefs) & !grepl("sq", coefs)) %>%
pull(value)),
"canopy_2" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("canopy_sq", coefs)) %>% pull(value)),
"slope" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("slope", coefs)) %>% pull(value)),
"herby" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("herby", coefs)) %>% pull(value)),
"memory" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("memory", coefs) & !grepl("sq", coefs)) %>%
pull(value)),
"memory_2" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("memory_sq", coefs)) %>% pull(value)),
"sl" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("step", coefs) & !grepl("log", coefs)) %>%
pull(value)),
"log_sl" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("log_step", coefs)) %>% pull(value)),
"cos_ta" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("turn_a", coefs)) %>% pull(value))
)
head(harmonics_df_2p)
## hour ndvi ndvi_2 canopy canopy_2 slope herby memory memory_2 sl log_sl
## 1 0.0 1.483291 -1.668939 0.4935256 -0.6770082 -0.2597942 0.06633518 2.172299 -1.254618 -0.8235112 0.05972889
## 2 0.1 1.492815 -1.675681 0.4900525 -0.6770977 -0.2628937 0.06622186 2.122835 -1.220851 -0.8577465 0.04886724
## 3 0.2 1.499561 -1.681229 0.4859317 -0.6766250 -0.2658112 0.06613404 2.072563 -1.186309 -0.8909046 0.03802145
## 4 0.3 1.503533 -1.685581 0.4811834 -0.6755962 -0.2685398 0.06607018 2.021594 -1.151077 -0.9229000 0.02721430
## 5 0.4 1.504745 -1.688739 0.4758296 -0.6740188 -0.2710729 0.06602851 1.970045 -1.115247 -0.9536502 0.01646880
## 6 0.5 1.503219 -1.690708 0.4698938 -0.6719016 -0.2734045 0.06600701 1.918035 -1.078908 -0.9830758 0.00580812
## cos_ta
## 1 -0.1083871
## 2 -0.1190373
## 3 -0.1295384
## 4 -0.1398633
## 5 -0.1499856
## 6 -0.1598791
# turning into a long data frame for plotting with ggplot
harmonics_df_2p_long <- pivot_longer(harmonics_df_2p, cols = !1, names_to = "coef")
Now we can plot all these together to see the temporal dynamics of the coefficients. Keep in mind that these are the scaled coefficients.
ggplot() +
geom_path(data = harmonics_df_2p_long,
aes(x = hour, y = value, colour = coef)) +
geom_hline(yintercept = 0, linetype = "dashed") +
scale_y_continuous(expression(Time-varying~parameter~values~beta)) +
scale_x_continuous("Hour", breaks = seq(0,24,2)) +
scale_color_discrete("Estimate") +
theme_classic() +
theme(legend.position = "bottom")
We can do the same with the natural scale coefficients.
harmonics_nat_df_2p <- data.frame(
"hour" = hour,
"ndvi" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("ndvi", coefs) & !grepl("sq", coefs)) %>%
pull(value_nat)),
"ndvi_2" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("ndvi_sq", coefs)) %>% pull(value_nat)),
"canopy" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("canopy", coefs) & !grepl("sq", coefs)) %>%
pull(value_nat)),
"canopy_2" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("canopy_sq", coefs)) %>% pull(value_nat)),
"slope" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("slope", coefs)) %>% pull(value_nat)),
"herby" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("herby", coefs)) %>% pull(value_nat)),
"memory" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("memory", coefs) & !grepl("sq", coefs)) %>%
pull(value_nat)),
"memory_2" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("memory_sq", coefs)) %>% pull(value_nat)),
"sl" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("step", coefs) & !grepl("log", coefs)) %>%
pull(value_nat)),
"log_sl" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("log_step", coefs)) %>% pull(value_nat)),
"cos_ta" = hour_harmonics_matrix %*%
matrix(coefs_clr %>% filter(grepl("turn_a", coefs)) %>% pull(value_nat))
)
head(harmonics_nat_df_2p)
## hour ndvi ndvi_2 canopy canopy_2 slope herby memory memory_2 sl log_sl
## 1 0.0 9.530534 -16.14346 2.099350 -4.041525 -0.2806850 0.1123308 24064517 -41415792435287 -0.002559160 0.04194061
## 2 0.1 9.570825 -16.20680 2.088156 -4.042115 -0.2843202 0.1121408 23612725 -40149347051813 -0.002672645 0.03863821
## 3 0.2 9.599658 -16.25926 2.074921 -4.039583 -0.2877461 0.1119922 23153804 -38855261640541 -0.002782553 0.03534126
## 4 0.3 9.617051 -16.30080 2.059711 -4.033965 -0.2909546 0.1118821 22688790 -37536794743147 -0.002888600 0.03205669
## 5 0.4 9.623059 -16.33146 2.042595 -4.025302 -0.2939379 0.1118074 22218746 -36197304753928 -0.002990512 0.02879152
## 6 0.5 9.617775 -16.35127 2.023650 -4.013644 -0.2966890 0.1117645 21744755 -34840240889741 -0.003088026 0.02555280
## cos_ta
## 1 -0.2087064
## 2 -0.2300259
## 3 -0.2510494
## 4 -0.2717232
## 5 -0.2919940
## 6 -0.3118097
# turning into a long data frame for plotting with ggplot
harmonics_nat_df_2p_long <- pivot_longer(harmonics_nat_df_2p, cols = !1, names_to = "coef")
Plot the natural scale coefficients. These are now on quite different scales as the covariates are on different scales - particularly the previous space use density, which we’ll omit from this plot.
ggplot() +
geom_path(data = harmonics_nat_df_2p_long %>%
filter(!coef %in% c("memory", "memory_2")),
aes(x = hour, y = value, colour = coef)) +
geom_hline(yintercept = 0, linetype = "dashed") +
scale_y_continuous(expression(Time-varying~parameter~values~beta)) +
scale_x_continuous("Hour", breaks = seq(0,24,2)) +
scale_color_discrete("Estimate") +
theme_classic() +
theme(legend.position = "bottom")
To update the Gamma and von Mises distribution from the tentative distributions from Fieberg et al. (2021), Appendix C, we just do the calculation at each time point.
tentative_shape <- 0.438167
tentative_scale <- 534.3507
tentative_kappa <- 0.1848126
hour_coefs_nat_df_2p <- harmonics_nat_df_2p %>% mutate(shape = tentative_shape + log_sl,
scale = 1/((1/tentative_scale) - sl),
kappa = tentative_kappa + cos_ta)
# turning into a long data frame
hour_coefs_nat_long_2p <- pivot_longer(hour_coefs_nat_df_2p, cols = !1, names_to = "coef")
Plot the updated movement parameters - notice that we’re scaling the scale parameter by 1000 to make it more interpretable.
ggplot() +
geom_path(data = hour_coefs_nat_long_2p %>%
filter(coef %in% c("shape", "kappa")),
aes(x = hour, y = value, colour = coef)) +
geom_path(data = hour_coefs_nat_long_2p %>%
filter(coef == "scale"),
aes(x = hour, y = value/1000, colour = coef)) +
geom_hline(yintercept = 0, linetype = "dashed") +
scale_y_continuous("Value of parameter") +
scale_x_continuous("Hour", breaks = seq(0,24,2)) +
scale_color_discrete("Estimate") +
theme_classic() +
theme(legend.position = "right")
Here we sample from the updated step length distribution (we could also follow the same process for turning angles) to generate a distribution for each hour of the day, and assess how well it matches the observed step lengths.
# summarise the observed step lengths
movement_summary <- buffalo_data %>% filter(y == 1) %>% group_by(id, hour) %>%
summarise(mean_sl = mean(sl), median_sl = median(sl))
## `summarise()` has grouped output by 'id'. You can override using the `.groups` argument.
# number of samples from the Gamma distribution
n <- 1e5
# create some empty variables to store the results
gamma_dist_list <- vector(mode = "list", length = nrow(hour_coefs_nat_df_2p))
gamma_mean <- c()
gamma_median <- c()
gamma_ratio <- c()
for(hour_no in 1:nrow(hour_coefs_nat_df_2p)) {
# sample from the Gamma distribution
gamma_dist_list[[hour_no]] <- rgamma(n, shape = hour_coefs_nat_df_2p$shape[hour_no],
scale = hour_coefs_nat_df_2p$scale[hour_no])
# summarise
gamma_mean[hour_no] <- mean(gamma_dist_list[[hour_no]])
gamma_median[hour_no] <- median(gamma_dist_list[[hour_no]])
gamma_ratio[hour_no] <- gamma_mean[hour_no] / gamma_median[hour_no]
}
gamma_df_2p <- data.frame(model = "2p",
hour = hour_coefs_nat_df_2p$hour,
mean = gamma_mean,
median = gamma_median,
ratio = gamma_ratio)
mean_sl_2p <- ggplot() +
geom_path(data = movement_summary,
aes(x = hour, y = mean_sl, colour = factor(id))) +
geom_path(data = gamma_df_2p, aes(x = hour, y = mean),
colour = "red", linetype = "dashed") +
scale_x_continuous("Hour", breaks = seq(0,24,2)) +
scale_y_continuous("Mean step length") +
scale_colour_viridis_d("Buffalo") +
ggtitle("Observed and modelled mean step length",
subtitle = "Two pairs of harmonics") +
theme_classic() +
theme(legend.position = "none")
mean_sl_2p
median_sl_2p <- ggplot() +
geom_path(data = movement_summary,
aes(x = hour, y = median_sl, colour = factor(id))) +
geom_path(data = gamma_df_2p, aes(x = hour, y = median),
colour = "red", linetype = "dashed") +
scale_x_continuous("Hour", breaks = seq(0,24,2)) +
scale_y_continuous("Median step length") +
scale_colour_viridis_d("Buffalo") +
ggtitle("Observed and modelled median step length",
subtitle = "Two pairs of harmonics") +
theme_classic() +
theme(legend.position = "none")
median_sl_2p
As we have both quadratic and harmonic terms in the model, we can reconstruct a ‘selection surface’ to visualise how the coefficient changes through time.
To illustrate, we have a coefficient for the linear term and a coefficient for the quadratic term for every hour of the day (or for every time point that we used to reconstruct the temporal dynamics, which is 0.1 hour in this case). Using these coefficients, we can plot the the quadratic selection curve at the scale of the environmental variable.
Using the natural scale coefficients of NDVI from the model, let’s take the coefficient from Hour 3 and plot the quadratic curve
# first get a sequence of NDVI values, starting from the minimum observed in the data to the maximum
ndvi_min <- min(buffalo_data$ndvi_temporal, na.rm = TRUE)
ndvi_max <- max(buffalo_data$ndvi_temporal, na.rm = TRUE)
ndvi_seq <- seq(ndvi_min, ndvi_max, by = 0.01)
# ndvi_seq
Similar to the harmonics, we add the linear and quadratic terms to get the quadratic curve
coef_hour <- 3
# we can separate to the linear term
ndvi_linear_coef <- hour_coefs_nat_df_2p$ndvi[which(hour_coefs_nat_df_2p$hour == coef_hour)]
ndvi_linear_coef
## [1] 6.748646
# multiply the coefficient by the values of NDVI
# the coefficient is positive so this will be a positive relationship
ndvi_linear_selection <- ndvi_linear_coef * ndvi_seq
plot(x = ndvi_seq, y = ndvi_linear_selection,
main = "Selection for NDVI - linear term at hour 3",
xlab = "NDVI", ylab = "Estimated selection")
lines(ndvi_seq, rep(0,length(ndvi_seq)), lty = "dashed")
# and the quadratic term - this term is negative so it will create a 'hat' rather
# than a 'bowl', with decreasing selection at the extremes
ndvi_quadratic_coef <- hour_coefs_nat_df_2p$ndvi_2[which(hour_coefs_nat_df_2p$hour == coef_hour)]
ndvi_quadratic_coef
## [1] -14.03231
# multiply the quadratic coefficient by the values of NDVI^2
ndvi_quadratic_selection <- ndvi_quadratic_coef * (ndvi_seq ^ 2)
plot(x = ndvi_seq, y = ndvi_quadratic_selection,
main = "Selection for NDVI - quadratic term at hour 3",
xlab = "NDVI", ylab = "Estimated selection")
lines(ndvi_seq, rep(0,length(ndvi_seq)), lty = "dashed")
# and sum them both together
ndvi_sum_selection <- ndvi_linear_selection + ndvi_quadratic_selection
plot(x = ndvi_seq, y = ndvi_sum_selection,
main = "Selection for NDVI - linear and quadratic terms at hour 3",
xlab = "NDVI", ylab = "Estimated selection")
lines(ndvi_seq, rep(0,length(ndvi_seq)), lty = "dashed")
We can see that the quadratic curve at Hour 3 shows highest selection for NDVI values slightly above 0.2, and the coefficient is greater than 0 for NDVI values from around 0 to 0.45.
For brevity we won’t plot the linear and quadratic terms separetely, but we can do so if needed.
Now for Hour 12
hour_no <- 12
# we can separate to the linear term
ndvi_linear_selection <- hour_coefs_nat_df_2p$ndvi[which(hour_coefs_nat_df_2p$hour
== hour_no)] * ndvi_seq
# plot(x = ndvi_seq, y = ndvi_linear_selection,
# main = "Selection for NDVI - linear term",
# xlab = "NDVI", ylab = "Estimated selection")
# and the quadratic term
ndvi_quadratic_selection <- (hour_coefs_nat_df_2p$ndvi_2[which(hour_coefs_nat_df_2p$hour
== hour_no)] * (ndvi_seq ^ 2))
# plot(x = ndvi_seq, y = ndvi_quadratic_selection,
# main = "Selection for NDVI - quadratic term",
# xlab = "NDVI", ylab = "Estimated selection")
# and the sum of both
ndvi_sum_selection <- ndvi_linear_selection + ndvi_quadratic_selection
plot(x = ndvi_seq, y = ndvi_sum_selection,
main = "Selection for NDVI - linear and quadratic terms at hour 12",
xlab = "NDVI", ylab = "Estimated selection")
lines(ndvi_seq, rep(0,length(ndvi_seq)), lty = "dashed")
Whereas for Hour 12, the coefficient shows highest selection for NDVI values around 0.5, and the coefficient is positive for all NDVI values above 0. The magnitude of selection is also greater than for Hour 3.
We can imagine viewing these plots for every hour of the day, where each hour has a different quadratic curve, but this would be a lot of plots. Another way to look at it is as a 3D surface, where the x-axis is the hour of the day, the y-axis is the NDVI value, and the z-axis (colour) is the coefficient value.
All we have to do is index over every time point and use the linear and quadratic coefficient values to reconstruct the quadratic curves, and then combine those curves into a 3D surface.
We use the same process as above, but now we index across every time point.
In the plotting we’ll add some contours for the 0 line, to show where the coefficient is positive and negative (conditional on all other covariates).
ndvi_min <- min(buffalo_data$ndvi_temporal, na.rm = TRUE)
ndvi_max <- max(buffalo_data$ndvi_temporal, na.rm = TRUE)
ndvi_seq <- seq(ndvi_min, ndvi_max, length.out = 75)
# Create empty data frame
ndvi_fresponse_df <- data.frame(matrix(ncol = nrow(hour_coefs_nat_df_2p),
nrow = length(ndvi_seq)))
# index across all time points
for(i in 1:nrow(hour_coefs_nat_df_2p)) {
# Extract the coefficient values for the linear and quadratic terms and
# multiply by the NDVI values
# Assign the vector as a column to the dataframe
ndvi_fresponse_df[,i] <- (hour_coefs_nat_df_2p$ndvi[i] * ndvi_seq) +
(hour_coefs_nat_df_2p$ndvi_2[i] * (ndvi_seq ^ 2))
}
ndvi_fresponse_df <- data.frame(ndvi_seq, ndvi_fresponse_df)
colnames(ndvi_fresponse_df) <- c("ndvi", hour)
ndvi_fresponse_long <- pivot_longer(ndvi_fresponse_df, cols = !1, names_to = "hour")
ndvi_contour_max <- max(ndvi_fresponse_long$value) # 0.7890195
ndvi_contour_min <- min(ndvi_fresponse_long$value) # -0.7945691
ndvi_contour_increment <- (ndvi_contour_max-ndvi_contour_min)/10
ndvi_quad_2p <- ggplot(data = ndvi_fresponse_long, aes(x = as.numeric(hour), y = ndvi)) +
geom_point(aes(colour = value)) + # colour = "white"
geom_contour(aes(z = value),
breaks = seq(ndvi_contour_increment, ndvi_contour_max, ndvi_contour_increment),
colour = "black", linewidth = 0.25, linetype = "dashed") +
geom_contour(aes(z = value),
breaks = seq(-ndvi_contour_increment, ndvi_contour_min, -ndvi_contour_increment),
colour = "red", linewidth = 0.25, linetype = "dashed") +
geom_contour(aes(z = value), breaks = 0, colour = "black", linewidth = 0.5) +
scale_x_continuous("Hour", breaks = seq(0,24,6)) +
scale_y_continuous("NDVI value", breaks = seq(-1, 1, 0.25)) +
scale_colour_viridis_c("Selection") +
# ggtitle("Normalised Difference Vegetation Index (NDVI)") +
theme_classic() +
theme(legend.position = "right")
ndvi_quad_2p
We can now see how the selection for NDVI changes throughout the day. We can clearly see the diurnal pattern of selection for NDVI, where the buffalo are seeking higher values of NDVI (centred on around 0.5) during the middle of the day, indicated by the positive coefficients, and relatively low selection for lower values of NDVI (centred on around 0.2) during the early morning and late afternoon, and negative coefficients for NDVI values above around 0.4, and strong selection against high values of NDVI. The selection against high NDVI values correlates with the movement dynamics, as they are likely moving through more open areas.
canopy_min <- min(buffalo_data$canopy_01, na.rm = TRUE)
canopy_max <- max(buffalo_data$canopy_01, na.rm = TRUE)
canopy_seq <- seq(canopy_min, canopy_max, length.out = 75)
# Create empty data frame
canopy_fresponse_df <- data.frame(matrix(ncol = nrow(hour_coefs_nat_df_2p),
nrow = length(canopy_seq)))
for(i in 1:nrow(hour_coefs_nat_df_2p)) {
# Assign the vector as a column to the dataframe
canopy_fresponse_df[,i] <- (hour_coefs_nat_df_2p$canopy[i] * canopy_seq) +
(hour_coefs_nat_df_2p$canopy_2[i] * (canopy_seq ^ 2))
}
canopy_fresponse_df <- data.frame(canopy_seq, canopy_fresponse_df)
colnames(canopy_fresponse_df) <- c("canopy", hour)
canopy_fresponse_long <- pivot_longer(canopy_fresponse_df, cols = !1, names_to = "hour")
canopy_contour_min <- min(canopy_fresponse_long$value) # 0
canopy_contour_max <- max(canopy_fresponse_long$value) # 2.181749
canopy_contour_increment <- (canopy_contour_max-canopy_contour_min)/10
canopy_quad_2p <- ggplot(data = canopy_fresponse_long, aes(x = as.numeric(hour), y = canopy)) +
geom_point(aes(colour = value)) +
geom_contour(aes(z = value),
breaks = seq(canopy_contour_increment, canopy_contour_max, canopy_contour_increment),
colour = "black", linewidth = 0.25, linetype = "dashed") +
geom_contour(aes(z = value),
breaks = seq(-canopy_contour_increment, canopy_contour_min, -canopy_contour_increment),
colour = "red", linewidth = 0.25, linetype = "dashed") +
geom_contour(aes(z = value), breaks = 0, colour = "black", linewidth = 0.5) +
scale_x_continuous("Hour", breaks = seq(0,24,6)) +
scale_y_continuous("Canopy cover", breaks = seq(0, 1, 0.25)) +
scale_colour_viridis_c("Selection") +
# ggtitle("Canopy Cover") +
theme_classic() +
theme(legend.position = "right")
canopy_quad_2p
The selection surface for Canopy cover tells a similar story to NDVI, with the buffalo selecting for higher values of canopy cover during the middle of the day, and lower values during the early morning and late afternoon, suggesting that they are seeking more cover during the hotter parts of the day. Overall the values of the coefficient are not as high as for NDVI, and only range between around -1 and +0.5, whereas for NDVI the coefficient ranges from -4 to +4, suggesting that NDVI is more influential.
memory_min <- min(buffalo_data$kde_memory_density, na.rm = TRUE)
# memory_max <- max(buffalo_data$kde_memory_density, na.rm = TRUE)
# take the 95% quantile instead of the max as the selection coefficients get quite large
memory_max <- quantile(buffalo_data |> filter(y == 1) |>
pull(kde_memory_density), probs = 0.95, na.rm = TRUE)
memory_seq <- seq(memory_min, memory_max, length.out = 75)
# Create empty data frame
memory_fresponse_df <- data.frame(matrix(ncol = nrow(hour_coefs_nat_df_2p),
nrow = length(memory_seq)))
for(i in 1:nrow(hour_coefs_nat_df_2p)) {
# Assign the vector as a column to the dataframe
memory_fresponse_df[,i] <- (hour_coefs_nat_df_2p$memory[i] * memory_seq) +
(hour_coefs_nat_df_2p$memory_2[i] * (memory_seq ^ 2))
}
memory_fresponse_df <- data.frame(memory_seq, memory_fresponse_df)
colnames(memory_fresponse_df) <- c("memory", hour)
memory_fresponse_long <- pivot_longer(memory_fresponse_df, cols = !1, names_to = "hour")
memory_contour_min <- min(memory_fresponse_long$value) # 0
memory_contour_max <- max(memory_fresponse_long$value) # 2.181749
memory_contour_increment <- (memory_contour_max-memory_contour_min)/10
memory_quad_2p <- ggplot(data = memory_fresponse_long, aes(x = as.numeric(hour), y = memory)) +
geom_point(aes(colour = value)) +
geom_contour(aes(z = value),
breaks = seq(memory_contour_increment, memory_contour_max,
memory_contour_increment),
colour = "black", linewidth = 0.25, linetype = "dashed") +
geom_contour(aes(z = value),
breaks = seq(-memory_contour_increment,
min(-memory_contour_increment, memory_contour_min),
-memory_contour_increment),
colour = "red", linewidth = 0.25, linetype = "dashed") +
geom_contour(aes(z = value), breaks = 0, colour = "black", linewidth = 0.5) +
scale_x_continuous("Hour", breaks = seq(0,24,6)) +
scale_y_continuous("Previous space use", labels = scientific) +
scale_colour_viridis_c("Selection") +
# ggtitle("Relationship to previous space use") +
theme_classic() +
theme(legend.position = "right")
memory_quad_2p
## Warning: `stat_contour()`: Zero contours were generated
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
## Warning: `stat_contour()`: Zero contours were generated
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
The relationship to previous space use is a bit more complicated to interpret, and for other more or less harmonics the trend changes, which is unlike NDVI and canopy cover, which are stable from 1 or more pairs of harmonics. There is consistent selection for higher previous space use density in the evening (between around 6pm and 10pm), suggesting that they are attracted to areas they have previously used, which may be foraging areas or similar, although this would be interesting to investigate in the field. It should also be noted that there is no 0 contour, as the entire selection surface is above 0, indicating that buffalo will, on average, always select for areas that have been previously. This selection for areas of previous space use lend the home ranging behaviour in the simulations. The coefficient values range from near 0 to up around +4, suggesting that previous space use is quite influential, particularly in the evening.
sessionInfo()
## R version 4.2.1 (2022-06-23 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_New Zealand.utf8 LC_CTYPE=English_New Zealand.utf8 LC_MONETARY=English_New Zealand.utf8
## [4] LC_NUMERIC=C LC_TIME=English_New Zealand.utf8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] cowplot_1.1.1 ggExtra_0.10.1 ggh4x_0.2.6 Rfast_2.0.7 RcppZiggurat_0.1.6
## [6] formatR_1.14 scales_1.2.1 glmmTMB_1.1.8 clogitL1_1.5 Rcpp_1.0.10
## [11] ecospat_3.5 TwoStepCLogit_1.2.5 survival_3.5-5 viridis_0.6.2 viridisLite_0.4.2
## [16] matrixStats_1.0.0 patchwork_1.1.2 ggpubr_0.6.0 adehabitatHR_0.4.21 adehabitatLT_0.3.27
## [21] CircStats_0.2-6 boot_1.3-28.1 MASS_7.3-59 adehabitatMA_0.3.16 ade4_1.7-22
## [26] sp_1.6-0 ks_1.14.0 beepr_1.3 tictoc_1.2 terra_1.7-23
## [31] amt_0.2.1.0 lubridate_1.9.2 forcats_1.0.0 stringr_1.5.0 dplyr_1.1.2
## [36] purrr_1.0.1 readr_2.1.4 tidyr_1.3.0 tibble_3.2.1 ggplot2_3.4.2
## [41] tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.3 tidyselect_1.2.0 lme4_1.1-32 htmlwidgets_1.6.2
## [5] grid_4.2.1 pROC_1.18.0 munsell_0.5.0 codetools_0.2-19
## [9] ragg_1.2.5 units_0.8-1 miniUI_0.1.1.1 withr_2.5.0
## [13] audio_0.1-10 colorspace_2.1-0 highr_0.10 knitr_1.42
## [17] rstudioapi_0.14 ggsignif_0.6.4 Rdpack_2.4 labeling_0.4.2
## [21] emmeans_1.8.5 TeachingDemos_2.12 bit64_4.0.5 farver_2.1.1
## [25] coda_0.19-4 vctrs_0.6.2 generics_0.1.3 TH.data_1.1-2
## [29] circular_0.4-95 xfun_0.39 timechange_0.2.0 randomForest_4.7-1.1
## [33] R6_2.5.1 isoband_0.2.7 cachem_1.0.7 reshape_0.8.9
## [37] promises_1.2.0.1 vroom_1.6.1 multcomp_1.4-23 nnet_7.3-18
## [41] gtable_0.3.3 mda_0.5-3 sandwich_3.0-2 rlang_1.1.0
## [45] systemfonts_1.0.4 splines_4.2.1 rstatix_0.7.2 TMB_1.9.10
## [49] earth_5.3.2 broom_1.0.4 checkmate_2.1.0 biomod2_4.2-2
## [53] yaml_2.3.7 reshape2_1.4.4 abind_1.4-5 backports_1.4.1
## [57] httpuv_1.6.9 Hmisc_5.0-1 tools_4.2.1 nabor_0.5.0
## [61] ellipsis_0.3.2 raster_3.6-20 jquerylib_0.1.4 RColorBrewer_1.1-3
## [65] proxy_0.4-27 plyr_1.8.8 base64enc_0.1-3 classInt_0.4-9
## [69] rpart_4.1.19 zoo_1.8-12 cluster_2.1.4 tinytex_0.48
## [73] magrittr_2.0.3 data.table_1.14.8 mvtnorm_1.1-3 fitdistrplus_1.1-8
## [77] mime_0.12 hms_1.1.3 evaluate_0.20 xtable_1.8-4
## [81] mclust_6.0.0 gridExtra_2.3 compiler_4.2.1 KernSmooth_2.23-20
## [85] crayon_1.5.2 minqa_1.2.5 htmltools_0.5.5 later_1.3.0
## [89] mgcv_1.8-42 tzdb_0.3.0 Formula_1.2-5 DBI_1.1.3
## [93] sf_1.0-12 Matrix_1.6-5 car_3.1-2 permute_0.9-7
## [97] cli_3.6.1 rbibutils_2.2.13 parallel_4.2.1 pkgconfig_2.0.3
## [101] numDeriv_2016.8-1.1 foreign_0.8-84 foreach_1.5.2 bslib_0.4.2
## [105] estimability_1.4.1 plotmo_3.6.2 digest_0.6.31 pracma_2.4.2
## [109] vegan_2.6-4 rmarkdown_2.21 htmlTable_2.4.1 PresenceAbsence_1.1.11
## [113] shiny_1.7.4 gtools_3.9.4 nloptr_2.0.3 lifecycle_1.0.3
## [117] nlme_3.1-162 jsonlite_1.8.4 carData_3.0-5 fansi_1.0.4
## [121] pillar_1.9.0 lattice_0.21-8 fastmap_1.1.1 plotrix_3.8-2
## [125] glue_1.6.2 gbm_2.1.8.1 iterators_1.0.14 bit_4.0.5
## [129] class_7.3-21 stringi_1.7.12 sass_0.4.5 maxnet_0.1.4
## [133] textshaping_0.3.6 poibin_1.5 e1071_1.7-13 ape_5.7-1